Particle-COMETS: From dynamic FBA to spatial simulations with particles¶

Eran Agmon, University of Connecticut

This notebook is dedicated to developing a particle-COMETS simulation step by step using Vivarium. Each section of this notebook incrementally adds layers to our model, illustrating the different computational methods and how they integrate.

In [1]:
import numpy as np
# spatio-flux's customized Vivarium
from spatio_flux import SpatioFluxVivarium, render_path

Make a Vivarium¶

view available types and processes

In [2]:
# make a fresh vivarium
vi = SpatioFluxVivarium()
In [3]:
# view the available types
vi.get_types()
Out[3]:
['tuple',
 'array',
 '',
 'length^1_5*mass^0_5/time',
 'length/time',
 'enum',
 'temperature',
 'mass^0_5/length^1_5',
 'length^2*mass/time^3',
 'mass/length^3',
 'length*mass/time^2',
 'mass/length*time^2',
 'length^2*mass/current^2*time^2',
 'current*time',
 'edge',
 'map',
 'interval',
 'printing_unit/length',
 'particle',
 'mass/time^2',
 'length*time/mass',
 'wires',
 'luminosity',
 'length^2/time^2',
 'emitter_mode',
 'current^2*time^4/length^2*mass',
 'current^2*time^3/length^2*mass',
 'length^3/mass*time^2',
 'mass/length*time',
 'length^3*mass/current^2*time^4',
 'substance/length^3',
 'composite',
 'union',
 'quote',
 'length^2*mass/current^2*time^3',
 'length',
 'length^1_5*mass^0_5/time^2',
 'length^4*mass/time^3',
 'current^2*time^4/length^3*mass',
 'length^2*mass/current*time^3',
 'time^2/length',
 'schema',
 'substrate_role',
 'current*time/substance',
 'length^2',
 'number',
 'length^2*mass/current*time^2',
 'float',
 'printing_unit',
 'length/time^2',
 'length^3/time',
 'boolean',
 'string',
 'length^2/time',
 'length^4*mass/current*time^3',
 'boundary_side',
 'length^2*mass/time',
 'mass^0_5/length^0_5*time',
 '/time',
 'bounds',
 '/temperature*time',
 'maybe',
 'tree',
 'length*temperature',
 'length^0_5*mass^0_5/time',
 'length^0_5*mass^0_5',
 'substance/time',
 'process',
 'mass/length',
 'protocol',
 'length^2*mass/time^2',
 'positive_float',
 'current*length^2',
 '/substance',
 'mass/temperature^4*time^3',
 'time/length',
 'time',
 'current*length*time',
 'current*time/mass',
 'list',
 'mass',
 'step',
 'reaction',
 'length^2*mass/substance*temperature*time^2',
 'length/mass',
 'integer',
 'length*mass/current*time^3',
 'path',
 '/printing_unit',
 '/length',
 'substance',
 'length^2*mass/temperature*time^2',
 'current',
 'kinetics',
 'current*time^2/length^2*mass',
 'any',
 'mass/time^3',
 'current*length^2*time',
 'luminosity/length^2',
 'length^3',
 'mass/current*time^2',
 'length*mass/current^2*time^2']
In [4]:
# view the available processes
vi.get_processes()
Out[4]:
['composite',
 'ram-emitter',
 'DynamicFBA',
 'json-emitter',
 'MinimalParticle',
 'console-emitter',
 'Particles',
 'DiffusionAdvection']

dFBA¶

Dynamic Flux Balance Analysis (dFBA) extends traditional Flux Balance Analysis (FBA) to model the dynamic behavior of metabolic networks over time, allowing for the simulation of growth and substrate consumption in a changing environment.

In [5]:
# inspect the DynamicFBA process config
vi.process_config('DynamicFBA')
Out[5]:
{'model_file': 'string',
 'kinetic_params': 'map[tuple[float,float]]',
 'substrate_update_reactions': 'map[string]',
 'biomass_identifier': 'string',
 'bounds': 'map[bounds]'}
In [8]:
# dfba_config  = vi.process_config('DynamicFBA', dataclass=True)  # TODO get dataclass to configure
dfba_config = {
    "model_file": "textbook",
    "kinetic_params": {
        "glucose": (0.5, 1),
        "acetate": (0.5, 2)},
    "substrate_update_reactions": {
        "glucose": "EX_glc__D_e",
        "acetate": "EX_ac_e"},
    "biomass_identifier": "biomass",
    "bounds": {
        "EX_o2_e": {"lower": -2, "upper": None},
        "ATPM": {"lower": 1, "upper": 1}}}
In [9]:
# make a fresh vivarium
v1 = SpatioFluxVivarium()

# add a dFBA process
v1.add_process(
    name="dFBA",
    process_id="DynamicFBA",
    config=dfba_config)
v1.diagram(dpi='70')
Out[9]:
No description has been provided for this image
In [10]:
mol_ids = ["glucose", "acetate", "biomass"]

# add the molecular fields
for mol_id in mol_ids:
    v1.add_object(
        name=mol_id,
        path=['fields'],
        value=np.random.rand())

v1.connect_process(
    name="dFBA",
    inputs={"substrates": {
                mol_id: ['fields', mol_id]
                for mol_id in mol_ids}},
    outputs={"substrates": {
                mol_id: ['fields', mol_id]
                for mol_id in mol_ids}})

# add an emitter to save results
v1.add_emitter()
v1.diagram(dpi='70', show_values=True)
Out[10]:
No description has been provided for this image
In [11]:
v1.set_value(path=['fields', 'glucose'], value=10)
v1.set_value(path=['fields', 'biomass'], value=0.1)
field = v1.get_value(['fields'])
print(field)
{'glucose': 10, 'acetate': 0.7464969100996626, 'biomass': 0.1}
In [12]:
# save a file with the exact simulation state
v1.save(filename='dFBA_t0')
Saved file: out/dFBA_t0.json
In [13]:
vx = SpatioFluxVivarium(document='out/dFBA_t0.json')
vx.run(interval=10)
vx.diagram(dpi='70', show_values=True)
Out[13]:
No description has been provided for this image
In [14]:
# run the simulation
v1.run(interval=60)
In [15]:
# view the timeseries
v1.get_timeseries(as_dataframe=True)
Out[15]:
/global_time /fields/glucose /fields/acetate /fields/biomass
0 0.0 10.000000 0.746497 0.100000
1 1.0 9.904762 0.758605 0.107943
2 2.0 9.802006 0.771581 0.116515
3 3.0 9.691146 0.785476 0.125766
4 4.0 9.571550 0.800341 0.135747
... ... ... ... ...
56 56.0 0.000000 0.000000 0.986275
57 57.0 0.000000 0.000000 0.986275
58 58.0 0.000000 0.000000 0.986275
59 59.0 0.000000 0.000000 0.986275
60 60.0 0.000000 0.000000 0.986275

61 rows × 4 columns

In [16]:
# plot the timeseries
v1.plot_timeseries(
    subplot_size=(8, 3),
    combined_vars=[
        [  # combine the variables into a single subplot
            '/fields/glucose',
            '/fields/acetate',
            '/fields/biomass'
        ]])
Out[16]:
No description has been provided for this image

Spatial dFBA¶

In [17]:
mol_ids = ["glucose", "acetate", "biomass"]
rows = 2
columns = 1

# make a fresh vivarium
v2 = SpatioFluxVivarium()
for mol_id in mol_ids:
    v2.add_object(
        name=mol_id,
        path=['fields'],
        value=np.random.rand(rows, columns))

# add a dynamic FBA process at every location
for i in range(rows):
    for j in range(columns):
        dfba_name = f"dFBA[{i},{j}]"
        v2.add_process(
            name=dfba_name,
            process_id="DynamicFBA",
            config=dfba_config)
        v2.connect_process(
            name=dfba_name,
            inputs={"substrates": {
                        mol_id: ['fields', mol_id, i, j]
                        for mol_id in mol_ids}},
            outputs={"substrates": {
                        mol_id: ['fields', mol_id, i, j]
                        for mol_id in mol_ids}})

# add an emitter to save results
v2.add_emitter()
v2.diagram(dpi='70')
Out[17]:
No description has been provided for this image
In [18]:
# change some initial values
v2.merge_value(path=['fields', 'glucose', 0, 0], value=10.0)
v2.merge_value(path=['fields', 'biomass', 0, 0], value=0.1)
field = v2.get_value(['fields'])
print(field)
{'glucose': array([[10.        ],
       [ 0.61243865]]), 'acetate': array([[0.32678828],
       [0.96336482]]), 'biomass': array([[0.1       ],
       [0.41108596]])}
In [19]:
# run a simulation
v2.run(60)
In [20]:
# get a list of all the paths so they can be plotted together in a single graph
all_paths = [
    [render_path(['fields', mol_id, i, j]) for mol_id in mol_ids]
    for i in range(rows)
    for j in range(columns)]

# plot the timeseries
v2.plot_timeseries(
    subplot_size=(8, 3),
    combined_vars=all_paths)
Out[20]:
No description has been provided for this image
In [21]:
v2.plot_snapshots()
No description has been provided for this image

Diffusion/Advection¶

This approach models the physical processes of diffusion and advection in two dimensions, providing a way to simulate how substances spread and are transported across a spatial domain, essential for understanding patterns of concentration over time and space.

In [22]:
bounds = (10.0, 10.0)
n_bins = (10, 10)
mol_ids = [
    'glucose',
    'acetate',
    'biomass'
]
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {'biomass': (0, -0.1)}

# make a fresh Vivarium
v3 = SpatioFluxVivarium()

# add fields for all the molecules
for mol_id in mol_ids:
    v3.add_object(
        name=mol_id,
        path=['fields'],
        value=np.random.rand(n_bins[0], n_bins[1]))

# add a spatial diffusion advection process
v3.add_process(
    name='diffusion_advection',
    process_id='DiffusionAdvection',
    config={
       'n_bins': n_bins,
       'bounds': bounds,
       'default_diffusion_rate': diffusion_rate,
       'default_diffusion_dt': diffusion_dt,
       'advection_coeffs': advection_coeffs},
    inputs={'fields': ['fields']},
    outputs={'fields': ['fields']})

# add an emitter to save results
v3.add_emitter()
v3.diagram(dpi='70')
Out[22]:
No description has been provided for this image
In [23]:
v3.run(60)
In [24]:
v3.show_video()

COMETS (Computation Of Microbial Ecosystems in Time and Space)¶

COMETS combines dynamic FBA with spatially resolved physical processes (like diffusion and advection) to simulate the growth, metabolism, and interaction of microbial communities within a structured two-dimensional environment, capturing both biological and physical complexities.

In [25]:
bounds = (20.0, 10.0)  # Bounds of the environment
n_bins = (12, 6)
mol_ids = ['glucose', 'acetate', 'biomass']
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {'biomass': (0, 0.1)}

# make a fresh vivarium
v4 = SpatioFluxVivarium()

# initialize the molecular fields
max_glc = 10
glc_field = np.random.rand(n_bins[0], n_bins[1]) * max_glc
acetate_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field[0:int(1*n_bins[0]/5), int(2*n_bins[1]/5):int(3*n_bins[1]/5)] = 0.1  # place some biomass

v4.add_object(
    name='glucose',
    path=['fields'],
    value=glc_field)
v4.add_object(
    name='biomass',
    path=['fields'],
    value=biomass_field)
v4.add_object(
    name='acetate',
    path=['fields'],
    value=acetate_field)

# add a diffusion/advection process
v4.add_process(
    name='diffusion_advection',
    process_id='DiffusionAdvection',
    config={
       'n_bins': n_bins,
       'bounds': bounds,
       'default_diffusion_rate': diffusion_rate,
       'default_diffusion_dt': diffusion_dt,
       'advection_coeffs': advection_coeffs    },
    inputs={'fields': ['fields']},
    outputs={'fields': ['fields']})

# add a dynamic FBA process at every location
for i in range(n_bins[0]):
    for j in range(n_bins[1]):
        dfba_name = f"dFBA[{i},{j}]"
        v4.add_process(
            name=dfba_name,
            process_id="DynamicFBA",
            config=dfba_config        )
        v4.connect_process(
            name=dfba_name,
            inputs={"substrates": {
                        mol_id: ['fields', mol_id, i, j]
                        for mol_id in mol_ids}            },
            outputs={"substrates": {
                        mol_id: ['fields', mol_id, i, j]
                        for mol_id in mol_ids}})

# add an emitter to save results
v4.add_emitter()
v4.diagram(dpi='70',
    remove_nodes=[f"/dFBA[{i},{j}]" for i in range(n_bins[0]-1) for j in range(n_bins[1])])
Out[25]:
No description has been provided for this image
In [26]:
v4.run(60)
In [27]:
v4.show_video()
In [28]:
v4.plot_timeseries(
    subplot_size=(8, 3),
    query=[
        '/fields/glucose/0/0',
        '/fields/acetate/0/0',
        '/fields/biomass/0/0'],
    combined_vars=[[
        '/fields/glucose/0/0',
        '/fields/acetate/0/0',
        '/fields/biomass/0/0']
    ])
Out[28]:
No description has been provided for this image

Particles¶

In [29]:
import numpy as np
from spatio_flux import SpatioFluxVivarium
In [30]:
vi.process_config('Particles')
Out[30]:
{'bounds': 'tuple[float,float]',
 'n_bins': 'tuple[integer,integer]',
 'diffusion_rate': {'_type': 'float', '_default': 0.1},
 'advection_rate': {'_type': 'tuple[float,float]', '_default': (0, 0)},
 'add_probability': 'float',
 'boundary_to_add': {'_type': 'list[boundary_side]',
  '_default': ['left', 'right']},
 'boundary_to_remove': {'_type': 'list[boundary_side]',
  '_default': ['left', 'right', 'top', 'bottom']}}
In [31]:
bounds = (10.0, 20.0)  # Bounds of the environment
n_bins = (20, 40)  # Number of bins in the x and y directions

v5 = SpatioFluxVivarium()

v5.add_process(
    name='particle_movement',
    process_id='Particles',
    config={
        'n_bins': n_bins,
        'bounds': bounds,
        'diffusion_rate': 0.1,
        'advection_rate': (0, -0.1),
        'add_probability': 0.3,
        'boundary_to_add': ['top']})
v5.connect_process(
    name='particle_movement',
    inputs={'fields': ['fields'],
            'particles': ['particles']},
    outputs={'fields': ['fields'],
             'particles': ['particles']})

v5.initialize_process(
    path='particle_movement',
    config={'n_particles': 2})

v5.add_emitter()
v5.diagram(dpi='70')
Out[31]:
No description has been provided for this image
In [32]:
v5.save('particle_movement')
Saved file: out/particle_movement.json
In [33]:
v5.run(100)
v5_results = v5.get_results()
In [34]:
v5.plot_particles_snapshots(skip_frames=3)
Saving GIF to species_distribution_with_particles.gif
In [35]:
v5.diagram(dpi='70')
Out[35]:
No description has been provided for this image
In [36]:
v5.save(filename='v5_post_run.json', outdir='out')
Saved file: out/v5_post_run.json

Minimal particle process¶

In [37]:
import numpy as np
from spatio_flux import SpatioFluxVivarium,render_path
from spatio_flux.processes.particles import get_minimal_particle_composition
In [38]:
bounds = (10.0, 20.0)  # Bounds of the environment
n_bins = (20, 40)  # Number of bins in the x and y directions

v6 = SpatioFluxVivarium()

# make two fields
v6.add_object(name='glucose',
              path=['fields'],
              value=np.ones((n_bins[0], n_bins[1])))
v6.add_object(
    name='acetate',
    path=['fields'],
    value=np.zeros((n_bins[0], n_bins[1])))

v6.add_process(
    name='particle_movement',
    process_id='Particles',
    config={
        'n_bins': n_bins,
        'bounds': bounds,
        'diffusion_rate': 0.1,
        'advection_rate': (0, -0.1),
        'add_probability': 0.3,
        'boundary_to_add': ['top']})
v6.connect_process(
    name='particle_movement',
    inputs={
        'fields': ['fields'],
        'particles': ['particles']},
    outputs={
        'fields': ['fields'],
        'particles': ['particles']})

# add a process into each particles schema
minimal_particle_config = {
    'reactions': {
        'grow': {
            'glucose': {
                'vmax': 0.01,
                'kcat': 0.01,
                'role': 'reactant'},
            'acetate': {
                'vmax': 0.001,
                'kcat': 0.001,
                'role': 'product'}
        }}}
particle_schema = get_minimal_particle_composition(v6.core, minimal_particle_config)
v6.merge_schema(path=['particles'], schema=particle_schema['particles'])

# add particles to the initial state
v6.initialize_process(
    path='particle_movement',
    config={'n_particles': 1})

v6.diagram(dpi='70')
Out[38]:
No description has been provided for this image
In [40]:
# v6.set_value('particles', {})
# v6.initialize_process(
#     path='particle_movement',
#     config={'n_particles': 1})

v6.run(60)
v6_results = v6.get_results()
In [41]:
v6.plot_particles_snapshots(skip_frames=3)
Saving GIF to species_distribution_with_particles.gif

add diffusion¶

In [52]:
document = v6.make_document()
v6x = SpatioFluxVivarium(document=document)
particle_id = list(v6x.composite.state['particles'].keys())[0]
# v6x.composite.composition['particles'] = {}
# v6x.composite.composition['particles']
# v6x.composite.state['particles'][particle_id]['minimal_particle']['_type'] = 'process'
v6x.diagram(dpi='70')
Out[52]:
No description has been provided for this image
In [50]:
document = v6.make_document()
v6x = SpatioFluxVivarium(document=document)
In [51]:
v6x.run(60)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[51], line 1
----> 1 v6x.run(60)

File ~/code/vivarium-interface/vivarium/vivarium.py:574, in Vivarium.run(self, interval)
    571 if not self.emitter_paths:
    572     self.add_emitter()
--> 574 self.composite.run(interval)

File ~/code/process-bigraph/process_bigraph/composite.py:800, in Composite.run(self, interval, force_complete)
    798 for path in self.process_paths:
    799     process = get_path(self.state, path)
--> 800     full_step = self.run_process(
    801         path,
    802         process,
    803         end_time,
    804         full_step,
    805         force_complete)
    807 # apply updates based on process times in self.front
    808 if full_step == math.inf:
    809     # no processes ran, jump to next process

File ~/code/process-bigraph/process_bigraph/composite.py:697, in Composite.run_process(self, path, process, end_time, full_step, force_complete)
    694 if state is None:
    695     breakpoint()
--> 697 update = self.process_update(
    698     path,
    699     process,
    700     state,
    701     process_interval
    702 )
    704 # update front, to be applied at its projected time
    705 self.front[path]['time'] = future

File ~/code/process-bigraph/process_bigraph/composite.py:636, in Composite.process_update(self, path, process, states, interval, ports_key)
    614 """Start generating a process's update.
    615 
    616 This function is similar to :py:meth:`_invoke_process` except in
   (...)
    631     ``store``.
    632 """
    634 states = strip_schema_keys(states)
--> 636 update = process['instance'].invoke(states, interval)
    638 def defer_project(update, args):
    639     schema, state, path = args

File ~/code/process-bigraph/process_bigraph/composite.py:303, in Process.invoke(self, state, interval)
    302 def invoke(self, state, interval):
--> 303     update = self.update(state, interval)
    304     sync = SyncUpdate(update)
    305     return sync

File ~/code/spatio-flux/spatio_flux/processes/particles.py:197, in Particles.update(self, state, interval)
    194 x, y = get_bin_position(new_position, self.config['n_bins'], self.env_size)
    196 # Update local environment values for each particle
--> 197 updated_particle['local'] = get_local_field_values(fields, column=x, row=y)
    199 # Apply exchanges and reset
    200 exchange = particle['exchange']

File ~/code/spatio-flux/spatio_flux/processes/particles.py:42, in get_local_field_values(fields, column, row)
     40 local_values = {}
     41 for mol_id, field in fields.items():
---> 42     local_values[mol_id] = field[column, row]
     44 return local_values

KeyError: (1, 30)
In [ ]:
v6x.plot_particles_snapshots(skip_frames=3)

Particle-COMETS¶

In [46]:
dfba_config = {
    "model_file": "textbook",
    "kinetic_params": {
        "glucose": (0.5, 1),
        "acetate": (0.5, 2)},
    "substrate_update_reactions": {
        "glucose": "EX_glc__D_e",
        "acetate": "EX_ac_e"},
    "biomass_identifier": "biomass",
    "bounds": {
        "EX_o2_e": {"lower": -2, "upper": None},
        "ATPM": {"lower": 1, "upper": 1}
    }}
In [47]:
bounds = (10.0, 20.0)  # Bounds of the environment
n_bins = (4, 8)  # Number of bins in the x and y directions

v7 = SpatioFluxVivarium()

v7.add_process(
    name='particle_movement',
    process_id='Particles',
    config={
        'n_bins': n_bins,
        'bounds': bounds,
        'diffusion_rate': 0.1,
        'advection_rate': (0, -0.1),
        'add_probability': 0.1,
        'boundary_to_add': ['top']})
v7.connect_process(
    name='particle_movement',
    inputs={'fields': ['fields'],
            'particles': ['particles']},
    outputs={'fields': ['fields'],
             'particles': ['particles']})

# add fields and diffusion process
mol_ids = ['glucose', 'acetate', 'biomass']
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {'biomass': (0, -0.1)}

# initialize the molecular fields
max_glc = 10
glc_field = np.random.rand(n_bins[0], n_bins[1]) * max_glc
acetate_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field[0:int(1*n_bins[0]/5), int(2*n_bins[1]/5):int(3*n_bins[1]/5)] = 0.1  # place some biomass

v7.add_object(name='glucose', path=['fields'], value=glc_field)
v7.add_object(name='biomass', path=['fields'], value=biomass_field)
v7.add_object(name='acetate', path=['fields'], value=acetate_field)

# add a spatial diffusion advection process
v7.add_process(
    name='diffusion_advection',
    process_id='DiffusionAdvection',
    config={
       'n_bins': n_bins,
       'bounds': bounds,
       'default_diffusion_rate': diffusion_rate,
       'default_diffusion_dt': diffusion_dt,
       'advection_coeffs': advection_coeffs},
    inputs={'fields': ['fields']},
    outputs={'fields': ['fields']})

# add a dynamic FBA process at every location
for i in range(n_bins[0]):
    for j in range(n_bins[1]):
        dfba_name = f"dFBA[{i},{j}]"
        v7.add_process(
            name=dfba_name,
            process_id="DynamicFBA",
            config=dfba_config)
        v7.connect_process(
            name=dfba_name,
            inputs={
                "substrates": {
                    mol_id: ['fields', mol_id, i, j]
                    for mol_id in mol_ids}},
            outputs={
                "substrates": {
                    mol_id: ['fields', mol_id, i, j]
                    for mol_id in mol_ids}})

v7.initialize_process(
    path='particle_movement',
    config={'n_particles': 2})

v7.add_emitter()

v7.diagram(dpi='70',
           remove_nodes=[
               f"/dFBA[{i},{j}]"
               for i in range(n_bins[0]-1)
               for j in range(n_bins[1])]
           )
Out[47]:
No description has been provided for this image
In [48]:
v7.run(60)
In [49]:
v7.plot_particles_snapshots(skip_frames=3)
Saving GIF to species_distribution_with_particles.gif